 
//   C++ version of mathematical basic routines

//   $Id: Math.hpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#ifndef _MATH_HPP
#define _MATH_HPP

// STL
#include <stdexcept>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
//#include <limits>
#include <cassert>
#include <cmath>

using namespace std;

class math_error : public logic_error
{
public:
  explicit math_error ( const string& p_arg ) : 
    logic_error ( p_arg ) { }
}; 

namespace MathConst {
  // some mathematical constants
  const double EE = 2.7182818284590452354;
  const double PI = 3.14159265358979323846;
  const unsigned HIGH_BIT = 0x80000000;
}

// template & inlines

inline bool even(int k) { return (k == 2*(k/2)); }

inline bool odd(int k) { return (k != 2*(k/2)); }

// BLAS routines
double dasum (int,double *);
int dscal (int,double,double *);
int idamax (int,double *);
int daxpy (int,double,double *,double *);

int PolyFit (double *,int,double *,double *,int,double *,char *);
int PolyReg ();
double RootBi ();	// roots by bisection method

// Numerical Recipes - modified routines
void balance (double **,int);
void elmhes (double **,int);
void hqr (double **,int,double *,double *);
int  polint (double *,double *,int,double,double *,double *);
int  spline(double *,double *,int,double,double,double *);
double splint(double *,double *,double *,int,double);

/* Eispack routines */
void balance (double **,int);		/* balance matrix */
void elmhes (double **,int);		/* transform to upper Hessenberg */
void hqr (double **,int,double *,double *);	/* QR algorithm for upper H */

#endif // _MATH_HPP
